function Gamma_sx = uig(s,x)

 
% Gamma(s)   = \int_0^{\infty} t^(s-1) e^(-t) dt ; gamma function
% gamma(s,x) = \int_0^{   x  } t^(s-1) e^(-t) dt ; lower incomplete gamma
% Gamma(s,x) = \int_x^{\infty} t^(s-1) e^(-t) dt ; upper incomplete gamma

 
% Gamma(s)   = Gamma(s,0)
% Gamma(s,x) = Gamma(s) - gamma(s,x)

 

 
Gamma_s      = gamma(s);

 
gamma_sx     = (Gamma_s)*gammainc(x,s);          % Matlab uses the "regularized" lower incomplete gamma

 
Gamma_sx     = Gamma_s - gamma_sx;